Identification of Signature Genes in the PD-1 Relative Gastric Cancer Using a Combined Analysis of Gene Expression and Methylation Data

Background The morbidity and mortality rates for gastric cancer (GC) rank second among all cancers, indicating the serious threat it poses to human health, as well as human life. This study aims to identify the pathways and genes as well as investigate the molecular mechanisms of tumor-related genes in gastric cancer (GC). Method We compared differentially expressed genes (DEGs) and differentially methylated genes (DMGs) in gastric cancer and normal tissue samples using The Cancer Genome Atlas (TCGA) data. The Kyoto Encyclopedia of Gene and Genome (KEGG) and the Gene Ontology (GO) enrichment analysis' pathway annotations were conducted on DMGs and DEGs using a clusterProfiler R package to identify the important functions, as well as the biological processes and pathways involved. The intersection of the two was chosen and defined as differentially methylated and expressed genes (DMEGs). For DMEGs, we used the principal component analysis (PCA) to differentiate gastric cancer from adjacent samples. The linear discriminant analysis method was applied to categorize the samples using DMEGs methylation data and DMEGs expression profiles data and was validated using the leave-one-out cross-validation (LOOCV) method. We plotted the ROC curve for the classification and calculated the AUC (area under the ROC curve) value for a more intuitive view of the classification effect. We also used the NetworkAnalyst 3.0 tool to analyze DMEGs, using DrugBank to acquire information on protein-drug interactions and generate a network map of gene-drug interactions. Results We identified a total of 971 DMGs in 188 PD-1 negative and 187 PD-1 positive gastric cancer samples obtained from TCGA. The KEGG and GO enrichment analysis showed the involvement of the regulation of ion transmembrane transport, collagen-containing extracellular matrix, cell-cell junction, and peptidase regulator activity. We simultaneously obtained 1,189 DEGs, out of which 986 were downregulated, while 203 were upregulated in tumors. The enriched analysis of the GO's and KEGG's pathways indicated that the most significant pathways included an intestinal immune network for IgA production, Staphylococcus aureus infection, cytokine-cytokine receptor interaction, and viral protein interaction with cytokine and cytokine receptor, which have previously been linked with gastric cancer. The compound DB01830 can bind well to the active site of the LCK protein and shows good stability, thus making it a potential inhibitor of the LCK protein. To observe the relationship between DMEGs' expression and prognosis, we observed 10 genes, among which were TRIM29, TSPAN8, EOMES, PPP1R16B, SELL, PCED1B, IYD, JPH1, CEACAM5, and RP11-44K6.2. Their high expressions were related to high risks. Besides, those genes were validated in different internal and external validation sets. Conclusion These results may provide potential molecular biological therapy for PD-1 negative gastric cancer.


Introduction
Gastric cancer (GC), also referred to as stomach cancer, occurs when malignant tumors grow within the inner lining of the stomach. It is the most common type of cancer in the world to afect the gastrointestinal system. Te morbidity and mortality rates for this form of cancer rank second among all cancers, indicating the serious threat it poses to human health, as well as human life [1]. Gastric cancer is also the second most diagnosed cancer in China, with approximately 403,000 new cases (281,000 males and 122,000 females) and 291,000 deaths in 2015. It is also the third leading cause of death in relation to all types of cancer [2,3]. Currently, commonly used therapies to treat gastric cancer include chemotherapy, radiotherapy, and surgery. However, many patients sufering from stomach cancer are often frst diagnosed when the disease has reached an advanced stage. Despite the available therapy options, the overall response to treatment remains dismal, with a 5-year survival rate of less than 30% [4]. As the molecular biology of this stomach cancer became clearer, immunotherapy was introduced as a new therapeutic strategy. Indeed, immune checkpoint inhibitors can improve the therapeutic efcacy of gastric cancer patients by activating the patients' immune systems and enhancing their antitumor immune response [5]. Te tissues afected by gastric cancer often contain several infltrating immune cells, such as natural killer (NK) cells or T and B cells. Te occurrences of these activated immune cells, along with the associated efector molecules at elevated levels, indicate a longer survival [6]. Besides, the overexpression of PD-L1 also indicates a dismal survival in people with gastric cancer [7]. However, we still lack effective biomarkers for the prognoses of patients with PD-1 negative gastric cancer. Terefore, the discovery and application of new biomarkers were crucial for improving the therapeutic efects and predicting the prognoses of patients with PD-1 negative gastric cancer.
Gastric cancer is considered a heterogeneous disease, which involves a multistep and multifactor process. Multiple cumulative epigenetic and genetic changes-such as tumor suppressor gene mutations and hypermethylation-turn normal cells into tumor cells. Tis results in tumorigenesis and the development of stomach cancer and afects the disease's biological behavior [8]. Over the past few years, research focusing on gastric cancer has increasingly been paying attention to the epigenetic mechanisms controlling this type of cancer, including long non-coding RNAs (noncoding ribonucleic acid), microRNAs, histone modifcation, and DNA methylation [9]. DNA methylation is a biochemical process and a considerably important epigenetic factor in tumorigenesis. It has been previously reported that methylation of CpG islands may induce tumorigenesis and aberrant methylation of CpG islands may afect the function of the tumor suppressor gene by changing the expression level of CpG islands [10,11]. Because of its stable and easyto-detect properties, DNA methylation can also be used as a good biomarker and may become a meaningful target for cancer diagnosis and treatment. Gastric cancer has shown high rates of DNA hypermethylation [12]. However, previous studies found that various aberrant methylation can lead to gene inactivation and gene silencing and promote the development of gastric cancer [13,14]. In Epstein-Barr virus (EBV) positive gastric cancer, methylation patterns of several tumor suppressor genes, including CHD1 and P16, have been altered and were considered essential tumorigenesis mechanisms [15]. Meanwhile, Helicobacter pylori infection and gene promoter hypermethylation was involved in multiple steps of carcinogenesis. Helicobacter pylori infection may result in methylation of the trefoil factor (TFF) family 2 and E-cadherin promoters. However, and despite extensive research, the clinical impact of these studies remains limited. Te exact mechanism by which DNA methylation induces gastric cancer remains unclear, and the drugs targeting these potential biomarkers were lacking.
We used numerous analyses of aberrant methylation and gene expression for diferential expression genes (DEGs) and diferential methylation genes (DMGs) [16][17][18]. We applied bioinformatics to recognize candidate genes and try to understand the genetic foundation of the disease. Terefore, we must use a combined approach when analyzing the gene expression profles chip and the gene methylation chip in gastric cancer. Using TCGA data, this study compared DNA methylation and DEGs in normal tissues and tissues suffering from gastric cancer. Trough enrichment analyses, we screened the determining genes and pathways in stomach cancer that might afect the development of this type of cancer. We also screened DrugBank to identify the potential drugs that targeted upregulated DMEGs. We identifed and confrmed 10 DMEG genes in diferent internal and external validation sets. Our study showed that gastric cancer cells contained 114 dysregulated DMEGs, which might be used as molecular biomarkers to aid in the early recognition of PD-1 negative gastric cancer. Tese results may provide potential molecular biological therapy for gastric cancer.

Data Source and Preprocessing.
We acquired the most recent stomach adenocarcinoma (STAD) expression profles, methylation data, and information on clinical followup through the TCGA GDC API on October 2, 2020. We also downloaded the expression profles data and survival information of the GSE30219 data set from GEO.
Te TCGA data were processed following these methods: (1) Only normal samples and primary tumor samples were retained (2) Based on the PD-1 gene expression level, tumor samples whose PD-1 level was lower than the mean of PD-1 level of normal samples were defned as PD-1 negative samples, while those whose PD-1 level was higher were considered PD-1 positive samples ( -UTR (3′-untranslated region), TSS1500, intergenic region, TSS200, TSS (transcriptional start site), 5′-UTR (5′-untranslated region), frst exon, and gene body), each gene's average β value was calculated according to the corresponding relationship. We then took the integrated gene methylation data to obtain diferentially methylated genes (DMGs) using the R limma package. Te methylation regions identifed as FDR <0.01. And delta β-values >0.3 were classifed as hypermethylated regions, while those identifed as FDR <0.01 and delta β-values <−0.3 were classifed as demethylated regions.

Expression Profle Data Analysis.
Te DEGs between tumor samples and normal samples were analyzed using the R limma package, and P values were adjusted using the BH method. TCGA data were log2 transformed, and DEGs were identifed as follows: genes with FDR < 0.01 and log2FC > 1 were considered upregulated genes, while genes with FDR <0.01 and log2FC <− 1 were considered downregulated genes.

Diferentially Expressed Methylated Genes in Diferent
Regions. To identify the association between expression profle and methylation, we calculated the intersection of DEGs and DMGs as diferentially methylated and expressed genes (DMEGs). We then divided them into four distinct groups, namely HyperDown, HypoDown, HypoUp, and HyperUp. Te detailed grouping criteria are shown in Table 1.

Functional Enrichment Analysis.
We conducted the Kyoto Encyclopedia of Genes and Genomes (KEGG), the gene ontology (GO) enrichment analysis, as well as the pathway functional annotations analysis on DMGs and DEGs using the clusterProfler R package to detect the crucial functions linked to the diferential genes and the essential biological processes and pathways involved.

Evaluation of the Markers of Methylation and Expression
Profle. We conducted the principal component analysis (PCA) to diferentiate gastric cancer and paracancer samples for DMEGs. We applied the linear discriminant analysis (LDA) when categorizing the samples using DMEGs expression profles data and methylation data, respectively, and the leave-one-out cross-validation (LOOCV) was applied to validate this analysis. We plotted the ROC curve for the classifcation and calculated the AUC value for a more intuitive classifcation efect.

Identifcation of Potential Target Terapeutic Drugs.
We identifed the drugs potentially targeting upregulated DMEGs by screening DrugBank using NetworkAnalyst 3.0, a web-based visualization platform that provides a comprehensive analysis and interpretation of system-level gene expression data. We used it to analyze DMEGs and used the DrugBank database to analyze protein-drug interactions to generate a network map of gene-drug interactions.
Based on drug-target pairs from DrugBank and the key PPI network in the string (with a threshold score of 400), we calculated the proximity of the drug to gastric cancer. Here, we defned S as the gastric cancer-associated gene set, DMEGs; D as the node degree of the gastric cancerassociated gene set in PPI; T as the drug target gene set; and d(s, t) as the shortest path of s node and t node (where s ∈ S was gastric cancer-related gene and t ∈ T was the drug target gene): Where ω was the weight of the target gene. If the target gene was one of a gene in the BPH-related gene set, the calculation method was ω � -ln (D + 1). Otherwise, it was ω � 0.
We generated a simulated reference distance distribution for the drug. To do so, we simply randomly selected a group of protein nodes in the network as the target of the drug. Te number of nodes was the same as the target size (R). Ten, we calculated the distance d (S, R) between the simulated drug targets (representing the simulated drug) and DMEGs and generated the simulated reference distributions after 10,000 random repetitions. Meanwhile, the mean and standard deviation of the reference distribution of μ d (S, R) and σ d (S, R) and the corresponding observed distance were converted into standardized scores, namely, the degree of closeness (Z): We found that whether it was DMEGs or our randomly selected gene set (Supplementary Table 1); when drug distances distribution was concentrated at 1 to 2, we performed multiple hypothesis tests based on the random data from reference and selected 26 drugs (Supplementary Table 2) with small distances and FDR <0.01, as a candidate drug set related to the DMEGs gene set.

PD-1 Expression and Immune Characteristics.
To observe the expression of PD-1 in human gastric cancer, we compared the PD-1 expression in normal versus tumor cells, as illustrated in Figure 1(a). Te distribution of PD-1 expression in tumor samples signifcantly difered from that in normal samples. Besides, we compared the microenvironment diferences between normal and tumor samples, the scores of immune cells' cytolytic activity (CYT) (Figure 1(b)), and seven diferent types of immune T cells (Figure 1(c)) obtained by ssGESA. We found out that normal tumor samples' CYT and T cell scores were lower than those of tumor samples, which revealed the immunosuppression in gastric cancer tumor samples.

Analysis of DMGs.
To identify diferential methylation in gastric cancer, we compared methylation data from 188 PD-1 negative tumor samples and 187 PD-1 positive ones. In this research, we focused on two gene bodies, TSS200 and TSS1500. Troughout the three regions, we identifed 971 diferentially methylated genes (FDR < 0.05, | delta β-values | > 0.3). In the volcano map, as shown in Figure 2(a), there were 183 hypermethylated genes and 347 demethylated genes in the gene body region, 98 hypermethylation and 185 demethylations in the TSS200 region, and 153 hypermethylation and 282 demethylations in the TSS1500 region. We found that the number of hypermethylation was about twice higher than that of demethylation in all three regions, as shown in Figure 2(b). We can also see that 26 hypermethylated genes occurred in all three regions, while 61 were found in two regions, and 234 were only found in one region, as shown in Figure 2(c).
Twenty demethylated hypermethylated genes were present in all three regions, while 124 were present in two, and 506 were only present in one, as depicted in Figure 2(d), which confrms that methylation was related to the region. To explore these DMGs' roles, we analyzed the GO functional enrichment and KEGG pathway. Tere were 335 biological processes, 24 cellular processes, and 19 molecular functions, as depicted in Figure 2(e), which illustrates the regulation of ion transmembrane transport, collagencontaining extracellular matrix, cell-cell junction, and peptidase regulator activity.

Analysis of DEGs.
We used the limma package to analyze the DEGs between 188 PD-1 negative and 187 PD-1 positive samples. We fnally obtained 1,189 diferentially expressed genes, out of which 986 were downregulated, while 203 were upregulated (FDR < 0.01, | FC | > 1.5). Te volcano map of the DEGs is depicted in Figure 3(a). Trough an unsupervised hierarchical cluster analysis of these DEGs, we found that diferentially expressed genes could distinguish tumor samples from normal ones, as illustrated in Figure 3(b). We performed the KEGG pathway and GO enrichment analyses of down-and upregulated genes using the clusterProfler R software package. Tese genes were enriched into 51 KEGG pathways, 54 cellular components, 754 biological processes, and 68 molecular functions. As illustrated in Figure 3(c), the most prominent of these pathways included cytokine-cytokine receptor interaction, S, viral protein interaction with cytokine and cytokine receptor, and intestinal immune network for IgA production, which is known to infuence the development of gastric cancer. Te result of the GO enrichment analysis of the DEGs is illustrated in Figures 3

Joint Analysis of DEGs and Methylated Genes.
We examined DEGs and DMGs in three diferent regions (gene body, TSS200, and TSS1500), and to further explore these genes' impact on stomach cancer, we determined the intersection of the DEGs and DMGs as DMEGs (diferentially methylated and expressed genes), which were believed to play a more determining role in promoting or inhibiting the development of gastric cancer. We obtained 78 DMEGs in the gene body region, as shown in Figure 4(a); 38 in the TSS200 region; and 52 in the TSS1500 region, as shown in Figures 4(b) and 4(c). Figure 4(d) shows the diferential methylation fold change and diferential expression fold change of these DMEGs. As we can see, we obtained fve genes with the greatest fold change of expression, among HyperDown 5

Analysis of DMEGs Genes.
We identifed 114 DMEGs and illustrated their distribution on chromosomes, as shown in Figure 5(a). Tere were 17 DMEGs on the chr1 chromosome and more than 8 on each chromosome of chr2, chr3, chr8, chr11, and chr12. Besides, we found that the methylation patterns of DMEGs were similar and consistent in similar gene regions. To explore the diferences in gene expression and DNA methylation patterns between tumor and normal samples, we constructed a linear discriminant classifcation model using the DMEGs gene expression profles and methylation data obtained from GeneBody, TSS200, and TSS1500 and carried out a PCA (principal component analysis) and ROC analysis (Figures 5(b) and 5(c), respectively). Te PCA's results showed that both gene expression profles and methylation data from diferent regions were able to classify the tumor and normal samples. We analyzed 114 DMEGs using the clusterProfler R software to conduct the KEGG pathway and the GO enrichment analysis and enriched 184 biological pathways. We selected the 10 most signifcant ones, displayed in Figure 5(d).

Potential Target Terapeutic Drugs.
We used the Net-workAnalyst 3.0 tool and DrugBank database to analyze the DMEGs gene for protein-drug interactions and found 10 genes that interacted with the drug, as shown in Table 2. At the same time, we thought that DMEGs were crucial genes related to gastric cancer, and drugs targeting these genes may have a greater impact on the occurrence of gastric cancer. Based on the description given in the methodology section,  Journal of Oncology we generated a simulated reference distance distribution for the drug. We found that whether DMEGs or randomly selected gene sets are used as samples, the drug distance is concentrated in 1 to 2. We performed multiple hypothesis tests based on the random data and selected the small distance and FDR. A total of 26 drugs with <0.01, which were used as the candidate drug set related to the DMEGs gene set, were obtained through analysis (Supplementary Figure 1). Based on the DMEGs obtained in step 6 and the candidate drug set, we selected the intersection and obtained the table shown above. We selected the AP-22408 with a signifcant distance of Lck to the gene set.
We took LCK-AP-22408 as an example of molecular docking analysis to clarify the binding model between the candidate drug and the target. We frst searched the PDB database and downloaded the LCK protein to the docking experiments and the 3D structure shown in Figure 6(a). With Autodock Vina, we analyzed the binding pattern of the target and the candidate drug. Te results showed that the compound DB01830 could bind well to the active site of the LCK protein with a docking score of -9.2 kcal/mol and produced favorable hydrogen bonding with the amino acid residues LEU251, GLU317, and MET319 of the LCK protein. It also had hydrophobic interaction with VAL259, ALA271, LEU371, VAL301, and ALA381, as well as S-Π interaction with GLU320 and GLU249, as shown in Figures 6(b) and 6(c). Tese important interactions suggested, to some extent, that the compound can closely bind to the LCK protein. In addition, Figure 6(d) showed the conformational changes of the compound DB01830 bound to the LCK protein during the 100 ns molecular dynamics simulation. Te conformations of this compound remained relatively stable and generally lower than 3Å, except for slight fuctuations in the frst 25 ns. Tese results indicate that the compound DB01830 could bind to the LCK protein stably. Terefore, this compound may be a potential inhibitor of the LCK protein.  Journal of Oncology training set, summarized the dimensionality reduction results, and counted the number of each probe appearing for each of the 1,000 times (Figure 7(a)). Ten genes can be observed for a combination of the maximum occurrence frequency (Table 3). Tese 10 genes, with diferent lambda variation coefcient trajectories, are depicted in Figure 7(b), and diferent lambda standard deviation distributions can be seen in Figure 7(c). Finally, the KM curve analysis showed that these genes could diferentiate the low-and high-risk groups. Finally, we obtained the following risk score formula: We calculated the RiskScore for each sample according to their expression levels and arranged the distribution of this RiskScore, as illustrated in Figure 7(d). Furthermore, we applied the R software package timeROC to analyze prognostic categorization efciency for 1 year, 3 years, and 5 years, as illustrated in Figure 7(e). Te area under the ROC curve (AUC) from the model was quite high, with the 5-year AUC above 0.7. We also performed a z-score for the RiskScore and used the R software to determine the cut-of value. We used Maxistat to group low-and high-risk samples and plotted the KM curve, as indicated in Figure 7(f ). Te results revealed a signifcant diference between the two samples with log rank P > 0.0001. Among those, 35 samples were categorized within the high-risk group, while 105 were classifed within the low-risk group. However, one of the samples was missing survival information.
To evaluate the 10-gene signature's predictive value, we applied the same models and constants to those used for the training set and verifcation set. Similarly, we calculated the RiskScore for each sample in relation to their expression level and arranged the RiskScore distribution, as illustrated in Figure 8(a). Moreover, we performed the analysis of the ROC of the prognostic categorization efciency for 1 year, 3 years, and 5 years, as illustrated in Figure 8(b). Te AUC of 3 years was above 0.87. We also conducted a z-score for RiskScore and determined the cut-of value using the R software package Maxistat to group the samples into a lowrisk group and a high-risk group and plotted the KM curve, as depicted in Figure 8(c). A signifcant diference was observed between the two groups with log rank P � 0.0067 and HR � 2.48. Among those, 30 samples were grouped into the high-risk group, while 15 of them were classifed into the low-risk group. However, two of the samples were missing survival information.
To evaluate the 10-gene signature's predictive value, we applied similar models and coefcients to those used in the training set to the PD-1 negative sample from TCGA.      2 2 2 2 2 2 2 2 2   4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4   6 6 6 6 6 6 6 6 6 6 6   8 8 8 8 8 8    Similarly, we calculated each sample's RiskScore according to their expression level and plotted the RiskScore distribution, as depicted in Figure 9(a). It is evident that the group with a high RiskScore had relatively shorter OS compared to that with a low RiskScore, indicating that high-RiskScore samples had a poorer prognosis. Moreover, we performed the analysis of the ROC of prognostic categorization efciency for 1 year, 3 years, and 5 years, as illustrated in Figure 9(b). Te 3-year AUC was 0.73, while that of 5 years was 0.67. We also conducted a z-score to evaluate the RiskScore, determined the cut-of value using the R software package Maxistat to group the samples into a low-risk group and a high-risk group, and plotted the KM curve, as depicted in Figure 9(c). We noticed a signifcant diference between the two groups with log rank P > 0.0001, while HR � 2.48. Among those, 66 samples were classifed into the high-risk group, while 119 were categorized into the low-risk group. However, three of the samples were missing survival information.

16
Journal of Oncology used the same methods to develop a prognostic risk model to assess patients' prognoses. Based on the results of the proportion of PD-1 negative obtained in outcome 1, we found that the expression of PD-1 in the normal samples was approximately equal to the quartiles of the expression of PD-1 in the gastric cancer samples. So we took the lowest one-quarter from the 433 samples of GSE84437, namely, the 108 samples with the lowest PD-1 expression. Each sample's RiskScore was calculated according to their expression level, and the RiskScore distribution was plotted, as depicted in Figure 10(a). It is evident that the group with a high RiskScore had relatively shorter OS than the one with a low RiskScore, indicating that high-RiskScore samples had a poorer prognosis. Te expression of nine distinct signature genes rose with the increase of the RiskScore. Furthermore, we used the R software timeROC to analyze the ROC of RiskScore prognostic categories. We performed the analysis of the ROC of prognostic categorization efciency for 1 year, 3 years, and 5 years, as illustrated in Figure 10(b). Te 5-year AUC was 0.69. We also conducted a zscore for the RiskScore and determined the cut-of value using the R software package Maxistat to group the samples into a low-risk group and a high-risk group. We also plotted the KM curve, as depicted in Figure 10(c). We observed a signifcant diference between the two groups with log rank P > 0.00011, while HR � 1.75. Among those, 40 samples were regrouped into the high-risk group, while 68 were classifed into the lowrisk group.

Discussion
Gastric cancer is the most common cancer that afects the gastrointestinal system and is ranked as the second leading   Journal of Oncology cause of cancer-related death in the world [19]. Multiple gene methylation is closely linked to the development and occurrence of gastric cancer. In this study, we combined and analyzed two diferent types of PD-1 negative gastric cancer gene chips using a bioinformatics analysis tool, to reveal the epigenetic and genetic mechanisms of gastric cancer. We found some hub genes, which provided new ideas for the diagnosis and treatment of gastric cancer.
We identifed 971 DMGs using 188 PD-1 negative tumor samples and 187 PD-1 positive samples of gastric cancer in TCGA. Te results of the KEGG and GO function enrichment analysis indicated that the DMGs were linked to biological processes, including the regulation of ion transmembrane transport, collagen-containing extracellular matrix, cell-cell junction, peptidase regulator activity, and so on. At the same time, we obtained 1,189 DEGs, out of which 986 were downregulated, while 203 were upregulated. Te GO enrichment and KEGG pathways examination indicated that the most signifcant top pathways included the cytokine-cytokine receptor interaction, Staphylococcus aureus infection, viral protein interaction with cytokine and cytokine receptor, and intestinal immune network for IgA production, which has been linked with gastric cancer. DMEGs were chosen from the intersection of the two sets, suggesting that such genes may play a greater role in promoting or inhibiting the development of gastric cancer. Some genes, such as POU2AF1 and IYD, can be seen in diferent regions. POU domain class 2-associating factor 1 (POU2AF1) was a known B-cell transcription coactivator. Tis gene was both expressed in lymphocyte cells and the whole-genome RNA sequencing of human airway epithelial cells. POU2AF1, as well as its related pathways, may be a therapeutic target for smoking-related airway diseases [20]. Tis gene promoted multiple myeloma through amplifcation or other mechanisms as an oncogene [21]. Te functional polymorphism of the POU2AF1 gene 3′-UTR was linked with an increased predisposition to lymphoma [22]. In addition, it has been known for being present in systemic autoimmune diseases, including multiple sclerosis and rheumatoid arthritis [23,24]. Iodotyrosin deiodinase (IYD) is an essential thyroid hormone enzyme for the iodination homeostasis invertebrates, which enables the efective synthesis of thyroid hormone [25,26]. We found a specifc mutation of IYD in patients with abnormal iodine metabolism and congenital hypothyroidism [27][28][29][30]. Currently, there is no tumor-related research focusing on these genes. As newly discovered gene targets, POU2AF1 and IYD might have a pertinent function in the development of stomach cancer and have yet to be further explored.
We have identifed 114 DMEGs and noticed the methylation patterns of these DMEGs were similar and consistent in close gene regions. We constructed a linear decision classifcation model based on DMEGs. PCA results showed that both gene expression profles and methylation data from diferent regions were able to separate normal and tumor samples with high accuracy. To clarify the binding model between the candidate drug and the target, we used LCK-AP-22408 as an example to study the docking analysis. Te results showed that the compound DB01830 could bind well to the active site of the LCK protein, which shows good stability and may be a potential inhibitor of the LCK protein.
TRIM29 was reported to have diferent efects on different types of tumors. For instance, its overexpression showed a poor prognosis in gastric cancer and the progression of a tumor, as well as bowel, pancreatic, bladder, lung, pancreatic, liver, thyroid, endometrial, and ovarian cancers [31][32][33][34][35][36][37][38]. However, TRIM29 appeared to have an inhibitory efect on other tumors, including breast and prostate cancers [39,40]. A meta-analysis has been conducted to clarify the predictive value of TRIM29 for diferent human malignant diseases. Te study eventually included 2,046 eligible patients, and the results suggested an important relationship between TRIM29's expression upregulation and a poor prognosis in patients with malignant tumors [40]. Overall, TRIM29 may have a pertinent function in the carcinogenesis of various human malignant tumors and be a useful biomarker for predicting the prognosis of patients.
TSPAN8 was a member of the TSPAN superfamily and played a crucial function in regulating leukocyte transport, wound repair, and angiogenesis in the physical-biological progress [41]. Over the last few years, an increasing amount of studies have shown that the overexpression of TSPAN8 is determining metastasis and the development of various tumors. [42][43][44][45] Te specifc targeting of monoclonal       antibodies in Span8, as an emerging cancer therapeutic target, is becoming increasingly important in cancer therapy [46]. Eomesodermin (EOMES) is a T-box transcription factor that promoted the diferentiation and function of cytotoxic lymphocytes, regulated natural killer (NK) cells mediating antiviral and antitumor responses, and participated in inducing the exhaustion of CD8+ T cells [47][48][49][50]. Eomes had diferent efects on diferent types of tumors. A high expression of Eomes was linked to short overall survival in patients with colorectal cancer [51]. However, it was also considered an independent good prognostic factor in patients with metastatic renal cell carcinoma [52]. Te level of Eomes methylation was also closely related to the tumor. EOMES showed tumor-specifc DNA hypermethylation in people with advanced bladder cancer [53,54]. In addition, the aberrant methylation of the Eomes gene promoter region resulted in its downregulation and hepatocellular carcinoma initiation and progression [55]. Tese studies showed that the occurrence and development of malignant tumors were closely related to the development of malignancies. In 133 cases of esophageal squamous cell carcinoma (ESCC), higher Eomes levels were associated with a poor clinical prognosis [56]. Eomes hypermethylation may provide an efective approach in the diagnosis of ESCC patients [57].
PCED1B-AS1 increased signifcantly in glioblastomas (GBM) tissues and was closely related to the tumor's grade and size. Te high PCED1B-AS1 survival time was shorter compared to that of the low PCED1B-AS1 group. Practical experiments illustrated that PCED1B-AS1's gene silencing repressed the proliferation of glioma cells and induced apoptosis, indicating that PCED1B-AS1 provided an auspicious biomarker for the prognosis, as well as drug targets of glioma [58,59]. It has been demonstrated that PCED1-B-AS1 can promote the proliferation, migration, and epithelial to mesenchymal transition (EMT) process, thus promoting the progression of clear renal cell carcinoma [60]. At the same time, PCED1B-AS1 can increase the function and expression of PD-L1 in hepatocellular carcinoma and induce the apoptosis of immunosuppressive cells [61].
Carcinoembryonic antigen-associated cell adhesion molecule 5 (CEACAM5), a highly glycosylated protein of the CEACAM family, increased its expression in the human breast, stomach, colorectal, and non-small-cell lung carcinoma cancers by promoting the proliferation and migration of cancer cells to promote the progression of the tumor [62][63][64][65]. CEACAM5 was recognized as a tumor biomarker and an indicator of cancer recurrence. CEA-CAM5 might become a potential target for a variety of cancer therapies.
Te results obtained in this study may be limited considering all the genes and pathways were based on the bioinformatics approach. As a result, this research lacks clinical samples that would have been useful to validate the data obtained. Moreover, we did not examine clinical gene expression profles and clinical gene methylation data in this study. Terefore, to enhance the reliability of our results, we need to conduct additional experiments, including a wide range of animal and cell experiments in future work.

Conclusions
In short, this study used the bioinformatics method to perform a combined analysis of PD-1 negative gastric cancer, using a gene methylation microarray and a gene expression microarray. Tis helped provide better insights into the molecular mechanism and pathogenesis of gastric cancer. We found out that the hub genes might represent molecular targets and diagnostic markers for accurate diagnosis and efective treatment of gastric cancer. Our study has also revealed a wealth of gastric cancer-specifc signatures, many of which may serve as drug targets and probable biomarkers for clinical applications in the future.

Data Availability
Te data used to support the fndings of this study can be obtained from the corresponding author upon reasonable request.

Conflicts of Interest
Te authors declare that the research was conducted in the absence of any commercial or fnancial relationships that could be construed as potential conficts of interest.

Authors' Contributions
Y.H. designed the study. L.E. performed data analysis, and L.S. wrote the manuscript. W.Z.G performed data collection, and G.F.F. supervised the manuscript. Te current manuscript has been read and approved by all authors.

Supplementary Materials
Supplementary Figure 1